package simulate;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;

import jdistlib.Poisson;



public class SingleVariateRegression {
	
	public static double ratio(float lambda,int a, int numRuns) {
		Poisson px = new Poisson(lambda*a);
		Poisson py = new Poisson(lambda);
		double ratio = 0;
		int count = 0;
		for (int i = 0; i < numRuns; i++) {
			double y = py.random();
			double x = px.random()/a;
			
			if (x > 0) {
				System.out.println(x+","+y);
				count +=1;
				ratio += y/x;
			}
		}
		
		return ratio/count;
	}
	
	
	
	
	
	
	/*** we assume half the events are generated by x ~ pois(a*lambda1), y ~ pois(lambda1).
	 *   and the other half similarily with lambda 2
	 * @param lambda1
	 * @param lambda2
	 * @param a
	 * @param numRuns
	 * @return
	 */
	public static double w(float lambda1, float lambda2, int a, int numRuns) {
		int perClusterRuns = numRuns/2;
		Poisson px1 = new Poisson(lambda1*a);
		Poisson px2 = new Poisson(lambda2*a);
		Poisson py1 = new Poisson(lambda1);
		Poisson py2 = new Poisson(lambda2);
		long sxy = 0;
		long sxx = 0;
		long sx = 0;
		long sy = 0;
		for (int i = 0 ; i < perClusterRuns; i++) {
			int x = (int) px1.random();
			int y = (int) py1.random();
			sxy += x * y;
			sxx += x*x;
			sx += x;
			sy += y;
		}
		
		for (int i = 0; i < perClusterRuns; i++) {
			int x = (int) px2.random();
			int y = (int) py2.random();
			sxy += x * y;
			sxx += x*x;
			sx += x;
			sy += y;
		}
		
		double meanx = sx/(double)numRuns;
		double meany = sy/(double)numRuns;
		double w = (sxy - numRuns*meanx*meany)/(sxx - numRuns*meanx*meanx);
		//System.out.println(sx+","+sy+","+w*a);
	
		return w;
	}
	
	public static double slope(Poisson px1, Poisson py1, Poisson px2, Poisson py2, int D, int n1, int n2) {
		int x2 = (int) px2.random();
		int x1 = (int) px1.random();
		int xy = 0;
		for (int i = 0; i < n1; i ++) {
			int y = (int) py1.random();
			xy += x1*y;
		}
		for (int i = 0; i < n2; i++) {
			int y = (int) py2.random();
			xy += x2*y;
		}
		
		double result = xy/(double)(n1*x1*x1+n2*x2*x2);
		return result;
	}
	
	/*** calculates the slope assuming we sample x and y from Poisson px1 and py1 n/2 times and then from px2 and py2 n/2 times ***/
	public static double slope4(Poisson px1, Poisson px2, Poisson py1, Poisson py2, int n) {
		double xy = 0;
		double xx = 0;
		for (int i = 0; i < n/2; i++) {
			int y = (int) py1.random();
			int x = (int) px1.random();
			xy += x*y;
			xx += x*x;
		}
		
		for (int i = n/2; i < n; i++) {
			int y = (int) py2.random();
			int x = (int) px2.random();
			xy += x*y;
			xx += x*x;
		}
		
		return xy/xx;
	}
	
	
	/*** calculates the slope assuming x's are the expectation of x for each crime rate and we sample n/2 y values. ***/
	public static double slope3(double x1, double x2, Poisson py1, Poisson py2, int n) {
		double xy = 0;
		for (int i = 0; i < n/2; i++) {
			int y = (int) py1.random();
			xy += x1*y;
		}
		
		for (int i = n/2; i < n; i++) {
			int y = (int) py2.random();
			xy += x2*y;
		}int y = (int) py1.random();
		xy += x1*y;
		double xx = n/2*x1*x1 + n/2*x2*x2;
		
		return xy/xx;
	}
	
	/*** calculates the slope assuming sampling and x value for each area and a y value for each area, day combination. ***/
	public static double slope2(Poisson px1, Poisson py1, Poisson px2, Poisson py2, int A, int D) {	
		int xy = 0;
		int xx = 0;
		
		for (int a = 0; a < A/2; a ++) {
			int x = (int) px1.random();
			for (int i = 0; i < D; i ++) {
				int y = (int) py1.random();
				xy += x*y;
				xx += x*x;
			}
			
		}
		
		for (int a = A/2; a < A; a ++) {
			int x = (int) px2.random();
			for (int i = 0; i < D; i ++) {
				int y = (int) py2.random();
				xy += x*y;
				xx += x*x;
				
			}
		}
		
		double result = xy/(double)xx;
		return result;
	}
	
	public static double mean(double[] vals) {
		double total = 0;
		for (double val: vals) {
			total += val;
		}
		total = total/vals.length;
		return total;
	}
	
	public static double std(double[] vals, double mean) {
		double variance = 0;
		for (double val: vals) {
			variance += (val - mean)*(val - mean);
		}
		variance = variance/vals.length;
		double std = Math.sqrt(variance);
		return std;
	}
	
	/*** run the simulation where we sample an x value for each area and a y value for each area,day combination. ***/
	public static void simulate2(){
		int totalExpected = 130;
		int D = 365;
		
		double [] lambda1List = {0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1};
		for (double lambda1: lambda1List) {
			double lambda2 = 2* lambda1;
			int A =(int) Math.round(totalExpected/(1.5*lambda1));
			
			Poisson py1 = new Poisson(lambda1);
			Poisson py2 = new Poisson(lambda2);
			Poisson px1 = new Poisson(lambda1*D);
			Poisson px2 = new Poisson(lambda2*D);
		
			int samples = 100;
			double[] wlist = new double[samples];
			
			for (int i = 0; i < samples; i ++) {
				double w = slope2(px1, py1, px2, py2, A, 365);
				wlist[i] = w*365;
			}
			
			double mean = mean(wlist);
			double std = std(wlist,mean);
			double stderr = std/Math.sqrt(wlist.length);
			System.out.println(lambda1+","+mean+","+stderr);
		}
	}
	
	
	public static void main(String[] args) throws IOException {
		int D = 365;
		double [] lambda1List = {0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1};
		for (double lambda1: lambda1List) {
			double lambda2 = 2* lambda1;
			
			Poisson py1 = new Poisson(lambda1);
			Poisson py2 = new Poisson(lambda2);
			Poisson px1 = new Poisson(lambda1*D);
			Poisson px2 = new Poisson(lambda2*D);
			
			int samples = 100;
			double[] wlist = new double[samples];
			
			for (int i = 0; i < samples; i ++) {
				
				double w = slope4(px1, px2, py1, py2, 1000000);
				wlist[i] = w*365;
			}
			
			double mean = mean(wlist);
			double std = std(wlist,mean);
			double stderr = std/Math.sqrt(wlist.length);
			System.out.println(lambda1+","+mean+","+stderr);
		}
		
		
		
		
		
	}
	
	
	
	
	public static double runOnce(float[] baseRates, int D, int numCells, int days) {
		double xx = 0;
		int xy = 0;
		for (int i = 0; i < baseRates.length; i++) {
			Poisson py = new Poisson(baseRates[i]);
			Poisson px = new Poisson(D*baseRates[i]);
			
			for (int s = 0; s < numCells; s ++) {
				double x =  px.random()/D;
				for (int d = 0; d < days; d ++) {
					
					int y = (int)py.random();
					xx += x*x;
					xy += x*y;
				}	
			}
		}
		double w = xy/xx;
		return w;
	}
	
//	public static void main(String[] args) throws IOException {
//		BufferedWriter out = new BufferedWriter(new FileWriter("/home/finn/phd/data/20140203/10_houses.txt"));
//		int numHouses = 100000;
//		float[] baseRates = {.000273973f,0.000136986f}; //per house base rate = 1/10years and 1/20years
//		int housesPerCell = 10; //approx 40 for a 200mx200m cell
//		float[] cellRates = {baseRates[0]*housesPerCell, baseRates[1]*housesPerCell};
//		int numCells = numHouses/housesPerCell;
//		int D = 365;
//		int days = 100;
//		System.out.println("num_cells:"+numCells);
//		
//		int N = 10000;
//		float wtotal = 0;
//		for (int i = 0; i < N; i++) {
//			double w = runOnce(cellRates, D, numCells,days);
//			out.write(String.valueOf(w));
//			out.newLine();
//			System.out.println(w);
//			wtotal +=w;
//		}
//		System.out.println(wtotal/N);
//		out.close();
//	}
	
	

}
